#include <fstream> 
#include <life/lifepoly/jacobi.hpp>
//#include <life/lifepoly/legendre.hpp>
#include <life/lifepoly/im.hpp>

#include <math.h>

using namespace Life;







int main(int argc, char** argv )
{

    const int N = 3;
    // Legendre polynomials of degree N
    Jacobi<N, double> p( 0.0, 0.0 );

    std::ofstream fichier;
    fichier.open("legendre.txt", std::ios::out); //on ouvrre le fichier en ecriture 

    for( double i=-1;i<=1;i=i+0.1) { fichier << i <<" "<<p(i)<<"\n"; }

    fichier.close();

    //---------------------------------------------------------------------//

    ublas::vector<double> wr(N), xr(N);
    details::gaussjacobi<N, double, ublas::vector<double>, ublas::vector<double> >( wr, xr);
  
//    std::cout<< "\n"<<JacobiBatchEvaluation<3,double>( 0., 0., xr )<<"\n\n";
    fichier.open("gaussPt.txt", std::ios::out);
    for (int i=0;i<N;++i) { fichier<<xr(i)<<" "<< wr(i)<<"\n"; }

    //---------------------------------------------------------------------//

    double integral = 0;
    for( int q = 0; q < N; ++ q) {
      integral += wr(q)*pow(xr(q),5) ;
    }
    std::cout << "integral=" << integral << "\n";

    //---------------------------------------------------------------------//

    Jacobi<0, double> p0( 0.0, 0.0 );
    Jacobi<1, double> p1( 0.0, 0.0 );
    Jacobi<2, double> p2( 0.0, 0.0 );

    IM<1,4> im;
    double integral00 = 0;
    double integral10 = 0;
    double integral20 = 0;
    double integral11 = 0;
    double integral12 = 0;
    double integral22 = 0;

    double normalisation0 = 0;
    double normalisation1 = 0;
    double normalisation2 = 0;
    for( uint q = 0; q < im.points().size2(); ++ q) { 
      normalisation0 += im.weight(q)*p0(im.points()(0,q))*p0(im.points()(0,q)) ;
      normalisation1 += im.weight(q)*p1(im.points()(0,q))*p1(im.points()(0,q)) ;
      normalisation2 += im.weight(q)*p2(im.points()(0,q))*p2(im.points()(0,q)) ;
    }
    normalisation0=math::sqrt(normalisation0);
    normalisation1=math::sqrt(normalisation1);
    normalisation2=math::sqrt(normalisation2);

    for( uint q = 0; q < im.points().size2(); ++ q) { 
      integral00 += im.weight(q)*p0(im.points()(0,q))*p0(im.points()(0,q)) ;
      integral10 += im.weight(q)*p1(im.points()(0,q))*p0(im.points()(0,q)) ;
      integral20 += im.weight(q)*p2(im.points()(0,q))*p0(im.points()(0,q)) ;
      integral11 += im.weight(q)*p1(im.points()(0,q))*p1(im.points()(0,q)) ;
      integral12 += im.weight(q)*p1(im.points()(0,q))*p2(im.points()(0,q)) ;
      integral22 += im.weight(q)*p2(im.points()(0,q))*p2(im.points()(0,q)) ;
    }

    std::cout<<" L0-L0 :" <<integral00/(normalisation0*normalisation0)<<"\n";
    std::cout<<" L1-L0 :" <<integral10/(normalisation1*normalisation0)<<"\n";
    std::cout<<" L2-L0 :" <<integral20/(normalisation2*normalisation0)<<"\n";
    std::cout<<" L1-L1 :" <<integral11/(normalisation1*normalisation1)<<"\n";
    std::cout<<" L1-L2 :" <<integral12/(normalisation1*normalisation2)<<"\n";
    std::cout<<" L2-L2 :" <<integral22/(normalisation2*normalisation2)<<"\n";


    //---------------------------------------------------------------------//
}

